#!/usr/bin/env python
# -*- coding: UTF-8 -*-

import math
ConstMeanRadiusOfEarthInKm=6372.797 
ConstNauticalMile = 1852.0 # Nautical mile in meters


def rad2deg(rad):
    """Convert radians to degrees. Return float radians."""
    return rad * 180 / math.pi 
def deg2rad(deg):
    """Convert degrees to radians. Return float degrees."""
    return deg * math.pi / 180  
def distance(lat1, lon1, lat2, lon2): 
    lat1 = deg2rad(lat1)
    lon1 = deg2rad(lon1)
    lat2 = deg2rad(lat2)
    lon2 = deg2rad(lon2)
    theta = lon1 - lon2
    dist = math.sin(lat1) * math.sin(lat2) \
           + math.cos(lat1) * math.cos(lat2) * math.cos(theta)
    dist = math.acos(dist)
    dist = rad2deg(dist)
    meters = dist * 60 * ConstNauticalMile
    return meters
 
def bearing(lat1, lon1, lat2, lon2):
    """Calculate bearing from lat1/lon1 to lat2/lon2. Return float bearing angle (0 <= bearing < 360)."""
    #NM gives error for some values
    
    lat1 = deg2rad(lat1)
    lon1 = deg2rad(lon1)
    lat2 = deg2rad(lat2)
    lon2 = deg2rad(lon2)
    try:
        bearingradians = math.atan2(math.asin(lon2 - lon1) * math.cos(lat2), \
                                    math.cos(lat1) * math.sin(lat2) \
                                    - math.sin(lat1) * math.cos(lat2) \
                                    * math.cos(lon2 - lon1))
        bearingdegrees = rad2deg(bearingradians)
        if (bearingdegrees < 0):
            bearingdegrees = 360 + bearingdegrees
    except:return -999
         
    return bearingdegrees


def anglediff(a, b):
    """Calculate the difference of two degree angles"""
    d = abs(a - b)
    if d > 180:
        d = abs(360 - d)
    return d 
def newlatlon(lat, lon, dist, bear):
    """Calculate new latitude and longitude from lat, lon, distance (meters) and bearing (degrees, 0 <= bearing < 360).Return (newlat, newlon)."""
    bearrad = deg2rad(bear)
    # Latitude difference in degrees
    latdiff = ((dist / ConstNauticalMile) * math.cos(bearrad)) / 60
    # Mean latitude
    midlat = lat + (latdiff / 2.0)
    # Longitude difference in degrees
    londiff = (dist / ConstNauticalMile * math.sin(bearrad) /
               math.cos(deg2rad(midlat))) / 60
    return (lat + latdiff, lon + londiff)


def estimatediff(lat1, lon1, speedkmph, secs, bear, lat2, lon2):
    # Not very good name for this function...
    """
    Estimate the 2nd point's coordinates using the 1st point's
    lat+lon, speed in km/h, seconds after and bearing.
    Return error in meters between estimated lat+lon and lat2+lon2.
    """
    dist = speedkmph / 3.6 * secs
    lat_est, lon_est = newlatlon(lat1, lon1, dist, bear)
    return distance(lat2, lon2, lat_est, lon_est) 

def DistanceFE(lat1, long1, lat2, long2, unit='Km'): # Distance on Flat Earth
    degrees_to_radians = math.pi/180.0 
    phi1 = (90.0 - lat1)*degrees_to_radians
    phi2 = (90.0 - lat2)*degrees_to_radians 
    theta1 = long1*degrees_to_radians
    theta2 = long2*degrees_to_radians 
    cos = (math.sin(phi1)*math.sin(phi2)*math.cos(theta1 - theta2) + math.cos(phi1)*math.cos(phi2))
    arc = math.acos( cos ) 
    if unit=='Km':return (arc * ConstMeanRadiusOfEarthInKm)
    else:return arc 

def geoDistFromTo  (fromLat, fromLng, toLat, toLng, lang, shortForm ):
    dst=DistanceFE(fromLat, fromLng, toLat, toLng)
    bear=bearing(fromLat, fromLng, toLat, toLng)  
    if bear !=-999 :card=degreesCardinal(bear, lang, shortForm)
    else: card=''
    return [dst, bear, card]
# -- Greek Stuff ----------------------------------------
def degreesCardinal(dg, lng='el', shortForm=False):
    if lng=='el':card=(u'Βόρεια', u'Βορειο-Ανατoλικά', u'Ανατολικά', u'Νοτιo-Ανατολικά', u'Νότια', u'Νοτιο-Δυτικά', u'Δυτικά', u'Βορειο-Δυτικά')
    dg=((dg+11.25)/45.0) 
    dg=int(dg)  % 8 
    tmp= card[dg]
    if shortForm:
        tmp=tmp.split('-')
        tmp=[s[0] for s  in tmp]
        return ''.join(tmp).encode('utf-8')
    else:   
        return tmp.encode('utf-8') 
def DistanceFEfromCenterOfGreece(lat1, long1):return DistanceFE(lat1, long1, 39.00, 22.00)
def geoDistfromAthensGreece  (lat1, long1, lang='el', shortForm=False):
    return geoDistFromTo(37.59, 23.44, lat1, long1, lang, shortForm )
